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ABSTRACT 


The  transient  behavior  of  the  phase  error  probability 
density  function  of  a  first-order  phase-locked  loop  in  the 
presence  of  additive  white  Gaussian  noise  is  determined 
using  numerical  techniques.   In  addition,  the  modulo-27T 
probability  density  function  and  the  cycle  slipping  density 
function  are  found.   The  algorithm  used  in  these  analyses 
involves  numerical  integration  of  a  normalized  and  factored 
Fokker-Planck  equation.   Results  are  shown  for  cases  in- 
volving various  signal-to-noise  ratios  and  initial  conditions. 
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I.   INTRODUCTION 

A.   BACKGROUND  INFORMATION 

The  space  age  spurred  great  interest  in  the  development 
of  phase-locking  techniques  since  fixed  tuned  receivers  were 
impractical  for  receiving  signals  experiencing  a  wide  range 
of  doppler  shifts.   Due  to  the  inherent  nonlinearity  of  the 
loop  phase  comparator,  early  attempts  to  analyze  phase- 
locked  loop  behavior  involved  the  use  of  graphical  phase 
plane  methods  which  are  summarized  by  Viterbi  [1].   The 
initial  use  of  Fokker-Planck  techniques  to  determine  the 
steady  state  probability  distribution  of  the  first  order 
loop  phase  error  was  accomplished  by  Tikhonov  [2],  [3].  Since 
then  much  work  has  been  done  using  the  Fokker-Planck  method 
of  analysis  but  little  has  been  published  dealing  with  the 
transient  behavior  of  the  phase-error  process. 

Various  techniques  have  been  employed  in  attempting  to 
statistically  describe  the  phase-error  process  of  first 
order  phase-locked  loops.   Stationary  phase-error  distribu- 
tions are  well  documented  by  Viterbi  [1]  and  Tikhonov  [2] 
and  [3].   Transient  solutions  of  the  Fokker-Planck  equation 
have  been  obtained  by  using  eigenfunction  expansions  [4] 
and  by  exploiting  the  similarities  between  the  Fokker- 
Planck  equation  and  a  one-dimensional  heat-flow  equation 
[5].   La  Frieda  [4]  arrived  at  a  statistical  description 
of  the  reduced  modulo-2-rT  phase-error  process  by  applying 


separation  of  variables.   This  application  yielded  an  eigen- 
value problem  which  he  solved  numerically.   Dominiak  and 
Pickholtz  ['5]  arrived  at  a  direct  numerical  solution  to 
the  transient  behavior  of  the  phase-locked  loop  (PLL)  by 
subjecting  the  Pokker-Planck  equation  to  the  same  numerical 
procedures  as  was  provided  for  the  one-dimensional  heat- 
flow  equation  by  von  Neumann  and  Richtmyer  [6].   Their 
resulting  partial  difference  equation  was  solved  using  a 
digital  computer. 

This  study  was  started  in  the  belief  that  the  previous 
work  in  the  analysis  of  the  transient  behavior  of  the  phase- 
locked  loop  could  be  expanded  and  improved.   La  Frieda's 
technique  is  applicable  only  for  the  reduced  modulo-2Tr 
solution  and  does  not  describe  the  actual  behavior  of  the 
phase-error  process.   Dominiak  and  Pickholtz  create  serious 
doubt  as  to  the  completeness  and  validity  of  their  technique 
in  that  their  conlcusion  of  insignificant  buildup  of  proba- 
bility mass  at  4>  =  ±k2TT  conflicts  greatly  with  qualitative 
estimates  of  Viterbi  [1]  and  Lindsey  [7]  and  with  one's 
intuition.   It  was  also  felt  that  Dominiak  and  Pickholtz 
did  not  make  a  good  choice  of  SNR  for  their  numerical  exam- 
ples.  The  stationary  probability  density  function  for  their 
examples  is  so  flat  that  any  buildup  of  probability  mass 
would  probably  be  very  small. 

In  this  paper  the  Pokker-Planck  equation  is  normalized 
using  several  constants  defined  for  the  purpose  [1]  and  [5]. 
The  resulting  equation  is  then  expressed  in  terms  of  the 


product  of  two  functions,  one  of  which  is  indepen- 
dent of  time.    The  other  function,   which  is  much 
smoother  than  the  original  probability  density,  is  found 
numerically  on  a  digital  computer  as  a  function  of  time. 
Multiplying  the  two  yields  the  solution  to  the  Fokker- 
Planck  equation  which  quantitatively  describes  the  transient 
behavior  of  the  phase-error  process  of  a  first  order  phase- 
locked  loop.   It  is  found  that  this  procedure  is  very  effi- 
cient computationally  and  readily  provides  numerical  results 
describing  the  phase  error  process,  the  modulo-2Tr  phase 
error  process,  and  the  cycle  slipping  statistics. 

B.   REVIEW  OF  PHASE-LOCKED  LOOPS 

A  phase-locked  loop  is  a  device  which  automatically 
controls  an  oscillator  or  periodic  function  generator  so  as 
to  operate  at  a  constant  phase  angle  relative  to  a  refer- 
ence signal  source.   Figure  1  shows  the  fundamental  hard- 
ware structure  of  the  loop  with  the  appropriate  signals. 
The  three  major  components  of  the  PLL  are: 

a)  A  multiplier  which  multiplies  the  input  signal  to 
the  loop  and  the  output  signal  of  the  voltage  controlled 
oscillator,  i.e. , 


x(t)  =  2AK1sine(t)cose ' (t) 


=  AK1{sin[e(t)-el(t)]  +  sin[e(t)  +  e'(t)]}, 


(1) 


where  e(t)  and  e!(t)  are  the  input  signal  phase  and  the 
voltage  controlled  oscillator  output  phase  respectively. 

b)  A  linear  low-pass  filter  with  zero-input  response 

e  (t)  and  impulse  response  f(t)  which  filters  out  the  sum- 
frequency  component  of  x(t)  while  passing  a  filtered  version 
of  the  difference-frequency  component,  i.e., 

t 
e(t)  =  e  (t)  +  /   x(t-X)f(x)dX   t  >_  0  .     (2) 
0      0 

c)  A  voltage  controlled  oscillator  (VCO),  the  output 
frequency  of  which  is  equal  to  the  sum  of  its  quiescent 
frequency  with  the  product  of  its  input  voltage  and  a 
proportionality  constant,  i.e., 

*$*!  =  Uq  +  K2e(t).  (3) 

Following  through  on  the  mathematics  [1]  it  can  be 

shown  that  for  e  (t)  =  0, 

o       ' 


d*' (t)  =w  +  K0  /  f(t-A)AKnsin[e(x)-e' (x)]dx. 
dt       o    d.   n         1 

0  (4) 

A 
For  phase  error  <fr(t)  =  e(t)-ef(t)  and  total  loop  gain 


df^tl  =  d^Ctj.  _  Wo  _  AK/  f(t-x)sin4>(x)dx.    (5) 
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Letting 


then 


and 


e1(t)  =  e(t)  -  toQt  (6) 

e2(t)  =  e'(t)  -  W()t,  (7) 


<j>(t)  =  e1(t)  -  e2(t)  (8) 


,  v   de  (t)        t 
2±^i-  =  — jfe —  -  AK  /  f(t-x)sin<{>(x)dx      (9) 


which  implies  the  PLL  model  shown  in  Fig.  2. 

Additive  noise  at  the  input  which  is  narrowband  Gaussian 

noise  n(t)  with  zero  mean  and  spectral  density  _o  near  the 

2 
signal  frequency,  can  be  expressed  as 


n(t)  =j2,{nn  (t)sinco  t  +  nn(t)cosu)  t}         (10) 

1         O       d  o 

where  n,(t)  and  ru(t)  are  independent,  zero  mean  Gaussian 
processes  of  identical  spectral  densities,  which  are  the 
same  as  the  spectral  density  of  n(t)  but  translated  downward 
in  frequency  so  as  to  be  centered  about  zero  frequency.   A 
derivation  similar  to  the  derivation  of  (9)  is  made  [1]  with 
the  noise  having  the  form 


n'(t)  =  -n  (t)sine2(t)  +  n2(t )cos62(t )      (11) 


where  n'(t)  may  be  treated  as  a  white  Gaussian  random 
process  with  spectral  density   o.   This  leads  to  the  model 
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of  the  PLL  in  Figure  3,  [1].   The  equation  of  operation  of 
the  phase-locked  loop  is  now  given  by 


.  ,  x    d9  (t)      t 

Q^tr;  =  — ±£ K/  f(t-X)[Asin<Kx)  +  n'U)]dX. 

°  (12) 
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FIGURE  1.   Fundamental  Hardware  Structure  of  the 
Phase-Locked  Loop. 
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FIGURE  3.   Phase-Locked  Loop  Modeled  With  Additive  Noise 
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II.   FOKKER-PLANCK  EQUATION 

Let  a  first-order  phase-locked  loop  have  an  Input  signal 
e^t)  of  constant  frequency  w,  i.e.,  e(t)=wt,and  a  voltage  controlled 
oscillator  with  quiescent  frequency  w  .   Additive  thermal 
noise,  n(t),  at  the  front  end  of  a  receiver  is  a  zero-mean 
Gaussian  process  whose  spectral  density  is  essentially  flat 
over  the  frequency  range  of  the  receiver  and  may  be  assumed 

to  be  white  with  value  _o.   The  equation  describing  the 

2 
loop  operation  [1]  is  from  (12) 


dtLti  =  (w_Wo)  _  K[Asin4>(t)  +  n'(t)]       (13) 


where  $(t)  is  "the  difference  between  the  incoming  signal 
phase  and  the  VCO  output  phase  and  K  is  the  total  loop 
gain. 

It  can  be  shown [1]  that#(t)  is  a  continuous  Markov  process 
The  noise  n'(t)  is  a  stationary  Gaussian  process  with  zero 
mean  and  it  is  nearly  white  over  the  frequency  range  of 
interest.   Therefore, 


N 
E[n'(t1)n'(t2)]  =  -f   6(t2-t1).  (15) 


Under  these  conditions,  the  probability  density  function 
p((|>,t)  of  the  Markov  process  phase  error  <i>(t)  must  satisfy 
the  Fokker-Planck  equation  [1],  i.e., 
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3*?U?t)    -  -  fj  C(«-«0-AKaint)p(#,t)] 


K2N 


g  3   P^>t}  (16) 


+  "5 

3<J> 


with  initial  condition 


p(4>,0)  =  «(♦-♦)  (17) 


and  normalizing  condition 

CO 

/  p(4>,t)d$  =  1  .  (18) 

It  is  convenient  to  express  (16)  in  normalized  form.   Let 

BL  =  AK/4      loop-noise  bandwidth  (19) 

2 

a  =  A  /N  BT    signal  to  noise  ratio  (20) 

o  L     & 

Y  =  (oj-u)  )/4BT   frequency  ratio  (21) 

O       Li 

t  =  *JB,t      normalized  time.  (22) 

Substituting  (19)  -  (22)  into  (16)  yields  a  normalized 
Fokker-Planck  equation  which  is 


jf^Sl  ■  |y  [(sln»-y)p(»,T)]  +  1  3  P^-T)    (23) 


An  analytic  solution  to  (23)  has  never  been  found. 

It  is  apparent  that  by  considering  only  very  small  phase 
errors,  i.e.,  sin<j>  *  $9  then  (13)  can  be  linearized.   Making 


lh 


the  same  approximation  in  (23)  gives  a  closed  form  solution 
which  satisfies  the  initial  and  boundary  conditions  on  p(<t>,x) 
This  solution  is 


pU»t)  = 


27r(l-e"2T) 


1/2 


exp 


a[<f>-Y-(4>0-Y)e   ] 
2(l-e-2T) 


which  is  Gaussian  with 


mean  =  m(-r)  =  y  +  (<j>   -  y)   e 


-T 


and 


2  /  -2t 

variance  =  o    (x)  =  (1  -  e   )/a 


(24) 


(25) 


(26) 


It  is  the  Gaussian  nature  of  {2k)   which  provides  the 
information  necessary  to  choose  the  initial  boundary  condi- 
tion function  for  the  numerical  algorithm  used  in  the  solu- 
tion of  (23).   However,  before  the  selection  of  an  initial 
condition  is  discussed,  it  is  appropriate  to  analyze  the 
technique  utilized  in  arriving  at  the  numerical  algorithm 
used  in  this  study. 
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III.   TRANSIENT  ANALYSIS 

A.   FUNDAMENTAL  TECHNIQUE 

The  fundamental  technique  employed  in  this  study  is  the 
treatment  of  the  normalized  solution  to  the  Fokker-Planck 
loop  equation  as  the  product  of  two  functions,  i.e., 

p(<J>,t)    =   gU,-r)    exp[a(cos<}>    +   y<J>)1    .  (27) 

This  exponential  factor  was  selected  on  both  mathematical 
and  intuitive  grounds.   Intuitively,  one  would  expect  that 
after  having  slipped  a  cycle  the  loop  dynamics  would  again 
exhibit  some  measure  of  steady  state  behavior  prior  to 
another  slip.   Lindsey  [7]  makes  a  similar  observation 
in  his  analysis  of  cycle-slipping  probabilities.   At  any 
one  time  the  loop  phase-error  tends  to  wobble  about  one  of 
an  infinite  set  of  points  spaced  by  2ir  radians.   It  seems 
reasonable  that  the  phase  error  probability  density  should 
resemble  a  mod-2Tr  stationary  density  function  in  each  2tt 
interval  with  a  broad  taper  away  from  the  origin  reflecting 
the  random  walk,  i.e.  cycle  slipping,  nature  of  the  problem. 

The  "envelope"  encompassing  this  succession  of  quasi- 
mod-2-n-  functions  can  be  easily  imagined  to  be  much  broader 
and  smoother  than  the  p(<}>,t)  function.  It  seems  plausible 
that  the  function  g(4),x)  from  (27)  is  this  smooth  "envelope". 
The  desirability  of  this  observation  is  that  the  "smoothness" 
of  g(4>,x)  would  allow  a  much  coarser  grid  to  be  used  in  a 
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numerical  integration  scheme  than  would  be  possible  with 
the  function  p  (<£,?)  . 

Applying  (27)  to  the  normalized  F-P  equation  (23)  one 
obtains 

;«(♦>*>  =  (Y-sin*)  ;g(*>T)  +  i  9V*'°  .      (28) 

3t  3<J>         a  9<j)2 

It  is  now  apparent  that  solving  for  g(4>,x)  is  computationally 
simpler  than  solving  for  p(<J>,t).   And,  since  the  exponential 
factor  is  a  constant  for  each  value  of  <}>  along  the  4>-axis, 
p(<J>,t)  can  be  reconstructed  readily  using  (27). 

Equation  (28)  is  evaluated  numerically  by  using  the 
central-difference  method  to  approximate  the  derivatives, 
such  that  with  A<j>  and  At  being  the  increments  in  <j>  and  t 
respectively,  then 

3g(4>,T)      a      g(j  ,T+AT)-g(4),T)  /2qv 

3t  At 


3g(^T)      a      g(«j)+A(j),T)-g((t.,T)  (30) 

34>  A<J» 


2 

9      gU,t)  rgU+Acfr  ,T)-g(4),T)  g(<|>  ,t)-g((j)-A4),T)-[      1_ 

3    2  L  A<fr  A<{>  J    A<j) 


„g(<fr  +  A<l>,T)    -    2gU,t)    +    gU-A<fr,T)  (31) 

(Acj))2 

and   it    follows   directly  that 

(32) 

g(*,T+Ax)=g(<J),T)+(Y-Sin<J))|~[g(<}>+A(J),T)-g(*-A<|),T)] 

+     ^— p      [g(*  +  A4),T)-2g(<(.,T)+g(<J>-A<J.,T)]. 

a     U«>r 
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B.   OTHER  CONSIDERATIONS 

1.   Initial  Condition 

It  is  the  response  of  (23)  to  an  input  impulse 
function  (17)  which  is  sought.   However,  it  was  felt  that 
a  finite  difference  method  algorithm  acting  on  the  extreme 
discontinuities  of  a  numerical  delta  function  of  a<j>  width  and 
1/A<j>  height  might  result  in  significant  errors  in  the  re- 
sults.  Therefore  the  Gaussian  nature  of  (24)  was  utilized 
in  choosing  an  initial  boundary  condition  for  the  algorithm. 
The  boundary  condition  selected  was  a  Gaussian  function 
which  corresponds  very  nearly  to  the  response  of  (23)  at 
some  very  small  time,  assuming  the  impulse  function  (17) 
existed  at  time  t  =  0.   To  check  the  accuracy  of  this  boun- 
dary condition  three  techniques  were  employed.   The  first 
was  to  plot  the  variance  of  the  calculated  phase  error 
versus  time  to  determine  whether  or  not  the  loop  is  opera- 
ting in  its  linear  region;  if  so,  the  variance  vs.  time 
curve  would  extrapolate  linearly  back  to  an  equivalent 
origin,  where  the  variance  is  zero  at  x  =  0 .   The  second 
was  to  actually  test  the  phase  error  response  to  the  numer- 
ical delta  function  mentioned  above.   The  third  was  to 
compare  stationary  and  transient  mod-27r  results  with  those 
already  published  [1],  [4]  and  [5], 

When  the  initial  condition  is  the  Gaussian  function,  the 
equivalent  time,  t  ,  at  which  this  condition  would  occur 
if  pU,0)  =  6(<j>-4>  )  is  given  by  (26)  as 
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„  2  =  „2(T  )  =  I  -   e"2T°  (33) 

O         O        a 


and  it  follows  that 


t   *  -  kln(l-aon2)].  (34) 

o      d.  o 


T 
O 


In  this  study  BT  t  =  0  and  B- 1  =  -jr-  are  the  times  at  which 
the  equivalent  delta  function  and  the  Gaussian  initial 
condition  existed  respectively.   In  [53 »  however,  BT t  =  0 
was  taken  to  be  the  time  of  existence  of  their  Gaussian 
initial  condition,  and  this  does  not  truly  represent  the 
dependency  of  p(<J>,t)  on  time. 
2.   End-value  Algorithm 

The  use  of  a  digital  computer  in  evaluating  this 
numerical  integration  scheme  requires  a  finite  interval 
which  comprises  the  sample  points  along  the  <j>-axis .   A 
potentially  large  source  of  error  in  any  computations  per- 
formed on  this  fixed  interval  is  the  method  by  which  the 
end  points  of  the  interval  are  evaluated  after  each  At  step 
The  central-difference  method  algorithm  in  this  study,  if 
applied  exclusively,  would  result  in  an  error  propagating 
inward  one  a<j>  increment  per  at  increment  starting  with  the 
first  at  increment.   This  error  would  make  the  procedure 
useless  for  large  t  unless  the  total  ^-interval  were  made 
prohibitively  large.   To  minimize  error  from  this  source, 
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a  sixth-order  Stirling  extrapolation  formula  is  used  to 
evaluate  both  of  the  end  points. 
3.   Stability 

Stability  in  a  numerical  algorithm  implies  that  the 
errors  generated  at  each  step  of  the  numerical  integration 
process  tend  to  decay  as  the  solution  progresses  in  time. 
The  method  used  in  this  study  to  determine  the  stability 
conditions  is  patterned  directly  after  von  Neumann's  Fourier 
series  technique  [6]  using  the  format  presented  in  [8].   The 
mesh  is  assumed  sufficiently  small  such  that  the  coeffi- 
cients in  (32)  may  be  treated  as  constants  over  a  region 
which  is  small  but  contains  many  mesh  points.   For  N  sample 
points  spaced  a<j>  radians  along  the  ^-interval,  the  initial 
errors,  i.e.  t  =  0,  which  are  assumed  to  be  located  at  the 
sample  points  can  be  expressed  for  each  point  in  complex 
exponential  form  as 


N-l    jB  k<A4>) 
Ek  =   E   Ane         »   k  =  0,1,2,. ,.,N-1    (35) 


where 


n=0 


6.  =  /„  ?w..v   .  (36) 


n  '  (N-1)(A<|>) 


Because  the  finite-difference  equations  will  always  be 
linear  and,  therefore,  separate  solutions  additive,  one  need 
only  consider  the  propagation  of  the  error  due  to  a  single 
term,  e.g.  e^   '  *  .   The  coefficient  A   is  a  constant  and 

'     to  n 
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can  be  neglected;  and  for  simplicity  the  subscript  n  on 
@  is  dropped. 

To  investigate  the  propagation  of  this  error  for 
increasing  x,  one  must  find  a  solution  of  the  finite-diff- 
erence equation  which  reduces  to  eJ       when  x  =  £  ( Ax )  =0 
r,   is  defined  as  the  error  at  the  k   point  along  the  <t> 
interval  at  the  l        iteration  along  x.   Assume  that  the 
error  has  the  form 


r         =    eJ6k(A<f>)    eXx    =      jBk(A<|>)       X£(At) 
kJl 

'  =  ei*H**)    ^  (37) 


where  £,   =   e      and  A,  in  general,  is  a  complex  constant. 
Clearly,  (37)  reduces  to  eJ'8k(A*)  when  l   =  0.   Also,  the 
error  will  not  increase  for  increasing  x  provided 


d  <  1.  (38) 


Since  the  error  r,  „  satisfies  the  same  difference 
equations  as  g(<f>,x),  substitution  of  (37)  into  (32)  and 
subsequent  division  by  e°  Y   £   yields 


C  =  1  +(Y-sin4>)^{exp[j8(A4»)]-exp[-j6(A$)]} 


+  -   At9  {exp[j6(A<}))]-2  +  exp[-jB(A*)]} 
a  (A*)2 

1+2  (Y-sin<J»)~rsin[ 6  (A<}> )]+(-) — A-^{cos[8(  A4>)  ]-l } 

U*)  (39) 
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A  T 

If  one  assumes  that  the  2(Y-sin<j>)— -sin[&  (A<f> )  ]  term  is 

Aq> 

negligible  for  large  N  and  very  small  A<|>,  then  from  (38) 
the  following  inequality  is  established: 


At  <  £(A4>)  .  (40) 


For  numerical   calculations   let  Ax  and  A<J>  be  related  as 
At  =  J(A<}>)2.  (4l) 

For  values  of  a  in  the  vicinity  of  unity  and  for  A4>  <<  1, 

then  At  <<  A<|>.   Hence,  the  dropping  of  the  term  containing 

At  .   .   ,  .  .„ .  , 

--—  is  justified. 
A4> 

4 .   Cycle  Slipping 

The  inclusion  of  absorbing  boundaries  around  some 
sub-interval  of  <j> ,  e.g.  [-tt,tt],  enables  one  to  determine  the 
probability  that  the  phase  error  has  not  exceeded  the  limits 
of  that  sub-interval  at  any  point  in  time.   As  the  density 
function  spreads  out,  any  probability  mass  which  reaches 
the  limits  of  the  sub-interval  is  removed  from  the  system 
entirely  by  forcing  the  density  function  to  be  zero  outside 
the  sub-interval.   The  effect  of  this  technique  is  to  make 
the  remaining  area  of  the  density  function  equal  to  the 
probability  that  the  phase  error  has  not  reached  the  limits 
of  that  sub-interval. 

The  cycle-slipping  density  function  is  denoted 
q($,T)  to  distinguish  it  from  the  corresponding  function 


-  22 


p(4>,x)  for  the  unbounded  case.   As  long  as  the  phase  error  of 

the  loop  remains  within  the  limits  of  a  specified  subinterval, 

the  density  function  qU,x),  with  initial  condition 

q(<f>,0)  =  <5(4>),  will  satisfy  the  Fokker-Planck  equation  (23), 

i.e. 

iaifex).  fT[(Sinj,-Y)q(4,,T)]  + 1 92q(t,o,     (il2) 


Following  the  development  from  Viterbi  [1],  for  y  =  $      =0 
the  probability  that  the  phase  error  <}>  has  not  reached  some 
value  <f>   by  time  t  is  given  by 


4> 
♦  (*)  =  ll   qU,x)d*  (H3) 


where 


qU£,0  =  q(-*A,T)  =  0   for  all  x  .        (44) 


The  fundamental  difference  between  q(<j>,x)  and  p  (<jj,t)  is 
that  q(<t>,x)  is  not  strictly  a  probability  density  function, 
unless  it  is  normalized  by  i|/(t).   The  numerical  algorithm 
used  to  generate  q(<J>,x)  is  the  same  as  for  generating  p(<j>,x) 
except  that  the  end  points  q(<t>  ,x)  and  q(-<|>  ,t)  are  forced 
to  be  zero.   The  end  points  for  p(<j>,x)  are  found  by  extrapo- 
lation, thus  simulating  an  infinite  interval. 

In  this  study  an  analysis  was  made  of  the  cycle 
slipping  probabilities  for  the  primary  cycle  slip,  i.e., 
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and  for  the  primary  half-cycle  slip,  i.e., 

[-♦t»#t]  ■  C-»/2»»/2]  .  (46) 

Plots  of  ^(t)  versus  x  show  the  time  dependent  behavior  of 
this  cycle  slipping  probability. 
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IV.   RESULTS 

The  entire  numerical  analysis  of  this  study  was  pro- 
grammed for  an  IBM-360/67  computer  using  Fortran  IV  language, 
level  G.   Double  precision  computation  was  used  exclusively 
to  minimize  round-off  and  truncation  errors.   The  initial 
condition  on  p(<j>,t)  in  (23)  was  a  Gaussian  function  with 

mean  4   and  standard  deviation  a   =  tt/20.   The  ^-interval  of 
yo  o  y 

investigation  was  [-llir,llir]  with  g(<f>sT)  and,  therefore, 
p(((),t)  being  computed  over  the  entire  interval  to  allow  for 
the  non-symmetrical  conditions  y  i-    0  and  <J>   ^0.   Stability 
considerations  specify  that  Ax  <  (a/2)(A<t>)  ,  and  to  accom- 
modate this  requirement  At  was  chosen  such  that 
At  =  (a/*0(A$)  .   The  sampling  interval  along  the  <}>-axis 
was  selected  to  be  A<t>  =  tt/50  except  in  the  cycle  slipping 
calculation  where  A<J>  =  tt/100.   For  the  examples  which  will 
be  shown  here,  it  was  found  that  negligible  improvements 

were  obtained  for  smaller  values  of  A<jj  and  a    .   In  fact, 

o  ' 

a  "delta  function"  of  width  A<}>  and  height  1/A<}>  was  used  as 
the  initial  condition  on  p(4>,t)  in  several  cases,  and  the 
results  obtained  varied  only  by  approximately  1  part  in 
20,000  from  those  obtained  with  the  initial  Gaussian  condi- 
tion above.   Table  I  indicates  the  various  conditions  which 
were  studied.   The  specific  values  were  chosen  so  as  to 
provide  comparison  with  previously  published  work  [1],  [^], 
[5]-   Also  included  in  the  table  are  the  accuracy 
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considerations  for  this  study  expressed  as  areas  of  the 
probability  density  function  after  the  indicated  number  of 
time  steps. 

TABLE  I 


a 

yo 

Y 

At 

Number  of 
Time  Steps 

Area  Under 
Curves 

0.1 

0 

0 

tt2x10"5 

20,000 

1.000 

1.0 

0 

0 

TT   XlO 

40,000 

1.000 

1.0 

0 

sin(Tr/4 ) 

TT   XlO 

20,000 

1.000 

2.0 

tt/2 

0 

2  -h 

2tt  xlO 

20,000 

1.000 

Figures  4,  5,  6,  and  7  are  graphs  relating  to  the  condi- 
tions a  =  1.0  and  y  =  4>0  =  0«   Figure  4  is  a  sample  of  the 
dynamics  of  the  g(<j>,x)  function  from  (27).   It  demonstrates 
the  "smoothness"  of  the  function  with  respect  to  the  p(<j>3x) 
function  which  is  shown  in  Figure  5.   Transitional  and 
stationary  modulo-2TT  densities  are  depicted  in  Figure  6. 
The  stationary  modulo-2Tr  density  function  is  seen  to  be  a 
reproduction  of  the  curve  determined  by  Viterbi  [1],   Of 
particular  interest  from  Figures  5  and  6  is  that  the  p(<f>,x) 
function  continues  to  spread  out  with  buildups  of  probability 
mass  at  the  ±k2TT  points  even  though  the  steady  state  modulo-2TT 
density  function  has  been  reached.   Figure  7  is  a  plot  of  the 
variance  of  the  phase  error  probability  density  function 
versus  time.   It  is  interesting  to  note  the  linear  nature  of 
the  curve  for  BT t  >  1.2. 
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La  Frieda's  eigenfunction  analysis  [4]  resulted  in 
transitional  and  steady  state  mod-2-rr  densities  which  were 
reproduced  using  the  algorithm  of  this  study.   Figure  8 
is  included  to  show  the  p(4>,x)  curves  corresponding  to 
each  of  the  times  used  in  [4]  with  an  additional  curve 
to  demonstrate  clearly  the  asymmetric  probability  mass 
buildup  associated  with  the  conditions  a  =  2,  y   =  0,  and 
4>   =  ir/2. 

Inasmuch  as  La  Frieda  used  a  different  normalizing 
scheme,  Dt,  in  his  study,  equivalent  times,  i.e.  BTt=DT/2, 
are  shown  on  the  graph  for  ease  of  reference.   The 
asymmetry  of  probability  mass  distribution  is  due  to  the 
main  body  of  p(i,x)  being  offset  from  <J>  =  0  for  some  finite 
period  of  time  while  it  is  shifting.   Since  p(cf>,x)  is 
settling  even  while  it  is  shifting,  more  mass  will  accumulate 
at  <j>  =  2tt  than  at  $  =  -2ir.   Once  the  main  body  of  p(<J>,x)  is 
settled  about  (J>  -  0,  however,  the  slippage  of  mass  in 
either  direction  becomes  equally  likely  and  the  level  of 
asymmetry  will  remain  constant  for  increasing  time. 

Viterbi  dealt  with  a  detuned  loop  which  was  also  ana- 
lyzed in  this  study.   For  this  case  a  =  1.0,  y  =  sin(ir/4), 
and  6=0  [1].   Figures  9  and  10  show  the  mod-2TT  prob- 
ability density  function  and  p(*,x)  for  these  conditions 
respectively.   As  before  the  steady  state  mod-2TT  density 
function  determined  by  this  study  is  a  reproduction  of  that 
found  by  Viterbi.   As  would  be  expected  and  is  demonstrated 
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in  these  two  curves,  the  detuned  loop  exhibits  a  definite 

tendency  to  slip  in  the  direction  of  the  steady-state  phase 

shift  caused  by  detuning. 

Figure  11  is  the  p(<j>,x)  function  corresponding  to 

a  =  0.1  and  v  ~  4>   =0.   This  set  of  conditions  is  one  of 
'     o 

those  selected  by  Dominiak  and  Pickholtz  [5]  in  their 
analysis  of  the  first  order  PLL  by  a  numerical  technique. 
This  curve  shows  the  buildup  of  probability  mass  which 
went  undetected  in  the  other  study.   It  is  felt  that  Domin- 
iak and  Pickholtz  did  not  find  this  buildup  for  several 
reasons.   Firstly,  they  did  not  look  look  enough  in  their 
time  frame.   Secondly,  their  attempt  to  deal  directly  with 
the  p(<J>,x)  function  rather  than  an  intermediate  function, 
e.g.  g(4>jO  in  this  study,  resulted  in  a  grid  mesh  size 
which,  although  mathematically  correct,  was  too  coarse  to 
detect  the  subtle  and  yet  significant  mass  slippage  of  the 
probability  density  function.   Finally,  from  [1]  it  is 
apparent  that  for  a  =  0.1  the  mod-2iT  density  function  is 
very  flat  and  is  not  a  good  choice  of  SNR  with  which  to 
search  for  probability  mass  buildup. 

Figures  12  and  13  show  the  dynamics  of  the  cycle  slip- 
ping density  functions,  q(<j>9x),  for  primary  cycle  slips, 
[-tt,tt],  and  primary  half-cycle  slips,  [-tt/2  ,  u/2]  .   These 
curves  are  plotted  for  the  conditions  a  =  0.5  and  <f>   =  y  =0 
In  Figures  1*1  and  15  are  the  curves  which  demonstrate  how 
the  cycle  slipping  probability,  iK  x ) ,  varies  with  time  for 
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the  cases  shown  in  Figures  13  and  14.   The  function  q(<j>,x) 
for  the  primary  cycle  slip  is  essentially  identical  to  the 
p(<t>jT)  function  until  a  significant  amount  of  probability 
mass  is  lost  at  the  absorbing  boundaries  at  <t>  =  ±ir.   Until 
this  loss  occurs  the  area  under  the  q(4>,x)  curve,  i.e.  ^(t), 
is  1.0.   The  time  at  which  i|>(t)  drops  below  some  value, 
e.g.  0.5,  can  be  determined  from  Figure  14.   From  Figures 
13  and  15  similar  data  may  be  obtained  for  the  primary 
half-cycle  slip  which  has  absorbing  boundaries  at  <J>  *  ±  —  . 
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FIGURE  A.   gU,x)  Density  Function  for  a  =  l .  0  ,  y  =  4'o=0  . 
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FIGURE  5.   Probability  Density  Function  of  Phase  Error, 
a=1.0,  y  =  4>o=0. 


31 


Lf\ 

in 

cu 

LTv 

vo 

vo 

rH 

CM 

CO 

O 

iH 

CM 

• 

• 

• 

O 

o 

H, 

/y 

II 

II 

a|                            > 

' 

-P 

-p 

•^                         .•/ 

J 

J 

hJ 

^          /v 

y 

PQ 

CQ 

CQ           /V 

y 

^-\__CTC 

6CC 

9GC 

\       tec 

I             zcr 

*■<*»£ 

\? 

' 

o 


'T3 

cti 

CO 


7  "e_ 

*\ 

o 

^  CD 

cd 
a, 


(qouT/sq.xun  z'o)    jojjs  asmjd  jo  jpd  -  (i*<|))d 
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FIGURE  7.   Probability  Density  Function  Variance  as  a 
Function  of  Time,  a=1.0. 


33 


o 


p 
ct 

p 
o 

52 

K 

re 

OJ       LT\ 

r-H       OJ       LO 

LTv 

• 

•H 

o     o     o 
II      II      II 

H        H        t- 
P     P     P 

OJ        J 
II          1 

P 

l\ 

VO      CO       rH 

oj    -=r     o 

vo     (\i     in 

O      H      OJ 

•                  •                 • 

o     o     o 
II      II      II 

p     p     p 
p-    p    p 

CQ      CQ      CQ 

CT\      ON 
OJ      CO 

•              • 

rH      CT\ 
II         II 
P      P 
P      J 

cq    m 

L         L 

^~~^^^^^Z 

i^£Er~- 

GTC 

(■z: 

9GC 

~^*-^C 

~^_    \ 
===^     ^^\ 

./ 

\ 

o 
C 

•H 
V. 

73 

o 


U 
o 
U 
U 
<D 

-,  <1> 

£  w 
"  ctJ 


(qoux/sq.xun  2*0)  J0O"J9  as^qd  jo  jpd  -  (i<i>)d 


FIGURE  8.   Probability  Density  Function  of  Phase  Error, 
a  =  2.0,Y  =  0,  4  =tt/2. 
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FIGURE  9.   Modulo-2^  Probability  Density  Function  of 
Phase  Error,  a=1.0,  y=sln(v/k) ,  $o=0* 
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FIGURE  10.   Probability  Density  Function  of  Phase  Error 
a  =  1.0,  Y  =  sin(TT/4),  4>o=0. 


36 


(qou-f/s^iun  z'Q)    jojj9  eseud  jo  jpd  -  (i'<t>)d 

FIGURE  11.   Probability  Density  Function  of  Phase  Error 
<x  =  0.1,  y  =  4>o=0. 
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FIGURE  12.   Primary  Cycle  Slip  Density  Function,  a=0.5. 
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FIGURE  14.   Primary  Cycle  Slipping  Probability  as  a 
Function  of  Time,  0=0.5. 
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FIGURE  15.   Primary  Half-Cycle  Slipping  Probability  as  a 
Function  of  Time,  a=0.5. 
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V.   CONCLUSIONS 

This  study  has  resulted  in  the  development  of  a  valid 
technique  for  numerically  finding  transient  solutions  for 
the  Fokker-Planck  equation  as  applied  to  the  first-order 
PLL.   That  the  technique  works  is  demonstrated  by  its  repro- 
duction of  existing  published  stationary  mod-2-rT  solutions 
in  [1]  and  transient  mod-27T  solutions  in  [4],   Results 
not  previously  available  for  the  phase-error  probability 
density  function  and  the  cycle  slipping  probability  are 
presented  here.   Most  important,  however,  is  that  the 
qualitative  notion  of  the  buildup  of  probability  mass 
around  points  spaced  2-rr  radians  has  been  verified. 
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